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1. Introduction 

In the present day, orbit determination by Global Positioning System (GPS) is not unusual. Especially for 
low-cost small satellites, position determination by an on-board GPS receiver provides a cheap, reliable and 
precise method. However, the original purpose of GPS is for ground users, so the transmissions from all of 
the GPS satellites are directed toward the Earth’s surface. Hence there are some restrictions for users above 
the GPS constellation to detect those signals. On the other hand, a desire for precise orbit determination for 
users in orbits higher than GPS constellation exists. For example, the next Japanese Very Long Baseline 
Interferometry (VLBI) mission “ASTRO-G” is trying to determine its orbit in an accuracy of a few 
centimetres at apogee. The use of GPS is essential for such ultra accurate orbit determination. This study 
aims to construct a method for precise orbit determination for such high orbit users, especially in High 
Elliptical Orbits (HEOs). There are several approaches for this objective. In this study, a hybrid method 
with GPS and an accelerometer is chosen. Basically, while the position cannot be determined by an 
on-board GPS receiver or other Range and Range Rate (RARR) method, all we can do to estimate the user 
satellite’s position is to propagate the orbit along with the force model, which is not perfectly correct. 
However if it has an accelerometer (ACC), the coefficients of the air drag and the solar radiation pressure 
applied to the user satellite can be updated and then the propagation along with the “updated” force model 
can improve the fitting accuracy of the user satellite’s orbit. In this study, it is assumed to use an 
accelerometer available in the present market. The effects by a bias error of an accelerometer will also be 
discussed in this paper. 


Table 1-1-1: Osculating Keplerian Orbital Elements 

EwTTu AA.nn.nn a i i\/i„,, oaa^ 


Epoch 

00:00:00.0, 16“ May, 2006 

Semi-major Axis 

19378.0 km 

Eccentricity 

0.62 

Inclination 

31.0 

RAAN 

0.0 

Argument of Perigee 

0.0 

Mean Anomaly 

0.0 



Figure 1-1-1: Image of ASTRO-G 

In this study, the next Japanese VLBI mission called "ASTRO-G" is used as a test mission. VLBI is a 
method of astronomical observation in the field of radio astronomy. The first VLBI mission was conducted 
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by Japanese satellite called "HALCA" which was the first VLBI satellite and finished its mission in 2005 
after 8 years operation. 

ASTRO-G is also planned by ISAS / JAXA, and it is scheduled to be launched in 2011. The orbit type of 
this mission is HEO, and the required orbit determination accuracy is a few centimetres at apogee. Figure 
1-1-1 shows the image of this mission. 

In Table 1-1-1 the osculating Keplerian orbital elements of ASTRO-G's orbit are summarised. As you can 
see in this table, the perigee altitude is 1,000 km while the apogee altitude is 25,000 km which is higher 
than the altitude of GPS orbits. The epoch time shown in Table 1-2-1 is the time set for this study. 

Figure 1-1-2 shows us the geometric relationship between the Earth, GPS satellite and ASTRO-G. 

GPS orbit 



Figure 1-1-2: Conceptual Diagram of GPS & ASTRO-G orbits 

In this study, it is intended to determine the orbit by using both of the GPS measurement data collected in 
the "lower" part of the orbit and the on-board ACC output data. The ACC allows to estimate more suitable 
force model, which will improve the orbit determination accuracy. 

2. Perturbation Forces 

In real situation, several perturbation forces are also applied to a satellite. When those forces are taken into 
account, the force vector a is described as; 

^ 3 ^ ® nonspherical ^3 -body(Sun) ^" 5 —body(Moori) ^ airdrag ^ SR 

r 

where, & nonsp i ienca i is the perturbation force cause by the non-spherical terms of the geo-potential model, 
a 3 -body(Sun) an d ^ 3-body(Moon > are f° rces by the third bodies, specifically the Sun and the Moon. 
a airdrag ' s ^ 1C force caused by the Earth's atmosphere and a SR is the solar radiation pressure. 

When we design the estimator utilizing ACC we have to consider the two disturbing forces i.e. the Air drag 
and the Solar radiation pressure. In the following sections those forces are introduced. 

2-1 Air Drag 

The disturbing force by the atmospheric drag has a strong influence against the motion of the LEO satellites. 
The next equation describes this perturbation force; 
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( 2 - 1 - 1 ) 
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where C D is the coefficient of drag which describe the susceptibility of the satellite to the air drag 
force. C D is approximately 2.2 in the upper atmosphere of the Earth when a flat plate model is applied. 
A is the cross-sectional area exposed to the atmosphere, and m is the mass of the satellite, p is the 
atmospheric density which differs along with the altitude. V rel is the velocity vector of the satellite 

relative to the atmosphere. This velocity vector is also described in the Earth Centred Inertial (ECI) 
coordinate system in this research. 


In this study, the following parameters are used; 

A = 50.o(m 2 ) m = 830.0(fcg) ( 2 - 1 - 2 ) 


These values are derived from HALCA's specification. 

Practically speaking, the atmospheric density p changes depending on the complicated interaction among 
some basic factors. However, in this research, simple values are used which changes with the altitude for 
simplicity. Figure 2-1-1 summarises the values of the mean atmospheric density [2]. And the derivative of 

“airdrag b y is Zer °- 

In this research the atmospheric density is interpolated linearly by using the values of Figure 2-1-1. The 
altitude above 1500km, it is defined that the density lineally reduces and reaches zero at the altitude 
1750km. 


For the simulation of this study the following description is used for simplicity. 

^ airdrag ^ drag \ ^ sat | Eva? 


( 2 - 1 - 3 ) 


The K drag appearing in the equation (2-1-3) is one of the state parameters estimated by the orbit 
estimator. 



2-2 Solar Radiation Pressure 

The radiation from the Sun also affects on the satellite motion as a perturbing force. The solar-radiation 
pressure is given by; _ . _ 

- _ PsR^R^Sun r sat-Sim 




m 

r 

| sat— Sun \ 
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( 2 - 2 - 1 ) 


where p SR is the solar pressure which equals to 4.5 1x10 6 N / m~ . C R is the reflectivity whose 
value is between 0.0 to 2.0. In this research the C R is set to 1.3. A Sun is the satellite’s exposed area to 
the Sun and m is satellite’s mass. These A Sun and m are same as the values shown in (2-1-2). 
r sat-Sun ' s l h e point vector from the satellite to the Sun. The derivative of a SR is zero. 


This perturbation force is also written in the following equation (2-4-2). The coefficient K SRP is 
estimated by the orbit estimator in this study. 


3. Visibility Analysis 


r 

Ty- sat- Sun 

U SR - ^ SRP I- 

\r 

| sat— Sun 


(2-4-2) 


3-1 Geometry 



Figure 3-1: Geometrical Allocation of GPS satellites and the user satellite 


To evaluate the feasibility of the orbit estimator designed in this study, the visibility analysis of GPS 
satellites should be investigated. Figure 3-1 shows the geometrical allocations of the GPS satellites and the 
user satellite. As can be seen, the user above the GPS orbits can only detect the signal from the GPS 
satellite which is beyond the Earth because all of the GPS signal transmission is directed toward the Earth, 
beside the beam width of the transmission antenna is not so wide. 


3-2 GPS Antenna Configuration 


To analyse the visibility of the GPS satellites in detail, the antenna patterns of the GPS transmission 
antenna and the GPS receiver's antenna are defined. As far as I know, the antenna configuration of the 
on-board GPS satellite antenna is not published. However some results of the investigation by GPS users 
are published. Basically, the antenna pattern for the GPS signal transmission should be designed to 
illuminate the Earth's surface from the altitude 20,200km. Hence the half angle should be about 13.88 
degrees. Figure 3-2-1 shows the antenna geometry. 


Orbit Altitude = 20,200 km 



Figure 3-2-1: The Geometry of the GPS transmission antenna 


13.88 deg 
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(a) GPS transmission antenna pattern 



Figure 3-2-2: Diagram of Antenna Patterns 


On the other hand, the beamwidth of GPS receiver's antenna is set to 150 degrees. Figures 3-2-2 shows the 
diagram of those antenna patterns. 


To figure out the visible number of GPS satellites, the shape of the satellite and the coordinate system was 
determined as shown in Figure 3-2-3. 

z 



Antenna Loading Facets 

Considered 6 antennas are loaded 
on all facets i.e. 

+ X & - X facets 
+ Y & - Y facets 
+ Z & - Z facets 


Figure 3-2-3: Definition of Coordinate System of the satellite 

As summarised in Figure 3-2-3, it is considered 6 GPS antennas are loaded on 6 panels of a cubed shape 
satellite. The satellite's body frame is fixed in the inertial space and the user satellite's shape is simply 
assumed a cube. Figure 3-2-4 shows the result of the visibility analysis. 
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Figure 3-2-4: Result of Visibility Analysis 

4. Orbit Estimator 
4-1 Kalman Filtering 

The orbit estimator designed in this study uses Kalman Filtering method. Kalman Filtering is a technique 
for computing the best estimate of the state parameters from imperfect observation values and uncertain 
dynamic model. Figure 4-1-1 shows the flow chart of this method. 
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In this study, the state parameters are the range and velocity vectors of the user satellite, and two 
coefficients of the air drag and the solar radiation pressure. The subscripts "drag" and "SRP" appearing in 
Figure 4-1-1 describe the air drag and the solar radiation pressure respectively. On the other hand, the 
observation values are the range, velocity and acceleration vectors. Note that the range and velocity vectors 
are described in ECI coordinate system while the acceleration vector's coordinate system is in the satellite's 
body frame. 



Figure 4-1-1: Flow Chart of Kalman Filtering 
4-2 Observation Matrix 


The observation matrix H is the matrix which describes the relationship between the state parameter X 
and the observation value Z . This matrix is defined by the following equation; 


In this study this matrix is; 

H _ az az az 
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(4-2-1) 


(4-2-2) 


Note that the acceleration a' appearing in the equation (4-2-2) is the force by the air drag and the solar 
radiation pressure only because the on-board ACC can only detect the non-conservative forces which are 
not proportional to the satellite's mass. Hence this force is described by; 


« = ~K d ,„ |5„ |v„ - K srp (4-2.31 

r 

| sat— sun | 

To distinguish this force from the acceleration of the force model the notation "dash" is used. By 
substituting the equation (4-2-3) into (4-2-2), then the following matrix is obtained; 
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where; 


1 3x3 @3x3 ®3x\ 

H — O ixi 1 3 x 3 0 3xl 

Ox.., M M , 


(4-2-4) 


= 


(4-2-5) 


M Jra g =-V sal V y 


(4-2-6) 


4-3 State Transition Matrix 


When we define the state transition matrix for 6 orbital parameters as cf* car , then this matrix can be 
obtained by solving the following equation; 


<J> = F<& 


(4-3-1) 


where; 


F _ a(r,v) _ 
d{r,v) 


dr dr 

dr dv 
dv dv 

dr dv. 


@3x3 1 3x3 

— o 

^3x3 

or 


(4-3-2) 


The vector a appearing in the equation (4-3-2) is the force vector shown in the section 3. The matrices 
0 3x 3 and f ]x3 are 3x3 zero matrix and the 3x3 identity matrix respectively. As can be seen in this 

equation we need to calculate the derivative of the force vector by the position vector of the user satellite to 
get the state transition matrix. 

When we describe the state transition matrix with four 3x3 matrices, the equation (4-3-1) is modified into 
the following equation; 

$J[0 3I 3 / 3 ,3 l.r*ii * 12 1 ( 4 - 3 -3) 

<*>21 <*>22 L D 0 3x3 \ l_<*> 2! <*>22. 


where D = (da/dr)l 3x3 . From the equation (4-3-3) we can obtain the state transition matrix C T> . On the 

other hand, the state transition matrix for the two coefficients of the air drag and the solar radiation 
pressure is simply described as following; 


O = 

^ ds 


(4-3-4) 


- 7 - 



because they are simply assumed constant. Consequently, the state transition matrix <!> in this study is 
written as; 

^car °6x2 (4-i 

o ct> 

\_^2x6 ^ ds J 


4-4 Process Noise Matrix 

The matrix Q is a process noise matrix which describes uncertainty of the dynamic model in the orbit 
estimator. According to [7], the formula of this matrix is provided as; 

h+ 1 h+ i r .. 

Q, = j (4 - 4 

h h 

where the vector it is a noise vector whose elements are white Gaussian noise. Here, we assume that the 
process noise emerges only in the vector W because the uncertainty of the model is in the force i.e. a . 
Hence the simplified orbit dynamics can be written; 

A _ 4 _ (4-d 

r = v, v = w 


then we obtain; 


r[u( b ! )n'(//)] = 


0,„ 


where S is the spectral amplitude associated with the white noise driving function, and S is Dirac 
delta function. U is specifically written as; 

“(#)=( 0 0 0 W M) w y(<5) (4 ~ 4 


S becomes " 1 " only when £ = T) . Therefore, the following equation is obtained from the equation 

(4-4-1); 
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4-5 Observation Noise Matrix 


The matrix R is an observation noise matrix which corresponds to the sensor's noise used in the observation. 
This matrix is described as following; 

M mg 0}x3 ®3x3 

R= 0„, M,„ 0 3 „ (4-5 

_ @3x3 @3x3 M acc _ 
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(4-5-2) 
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(4-5-3) 


(4-5-4) 


5. Simulation 

5-1 Simulation Configuration 

To figure out the contribution of an on-board ACC some simulations are carried out. Figure 5-1-1 shows 
the diagram of the simulation configuration. 



Figure 5-1-1: Diagram of Simulation Configuration 

As Figure 5-1-1 shows the reference orbit called "true orbit" is prepared at first. In real situation nobody 
knows the "true" orbit. The force model for the true orbit is summarised in Table 5-1-1. 

The geo-potential model used in this study is WGS 84 (EGM 96). Two disturbing forces by the Sun and the 
Moon are included. The coefficient of the air drag changes from 0.0 to 1.99e-13 because of the atmospheric 
density model applied to this study. On the other hand the coefficient of the solar radiation pressure is set to 
constant number. Note that the orbit estimator simply assumes that those coefficients of the air drag and the 
solar radiation pressure are constant. 

However, the real coefficients are not constant. They change depending on the solar activity whose 
magnitude is usually unpredictable, although the cycle of the period of the solar maximum and minimum is 
not unknown. 
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Table 5-1-1: Force Model of True Orbit 


Geo-potential Model 

36x36 

3-Body Disturbing Forces 

Sun & Moon 

Coefficient of Air Drag 

0.0 ~ 1.99e-13 

Coefficient of Solar Radiation Pressure 

3.53e-10 


5-2 Simulation Cases 

5-2-1 Casel: Influence of GPS Positioning Accuracy 

Here we can see the influence of the GPS positioning accuracy to the orbit determination error. Table 
5-2-1- 1 shows the relationship between the noise configurations of GPS positioning and the orbit 
determination errors when the orbit is determined only by the GPS measurement data. The initial 
coefficient of the air drag is set to 0.0 and the coefficient of the solar radiation pressure is 7.06e-10 which is 
twice larger than the true value. 


Table 5-2-1-1: GPS Positioning Accuracy and 

1 Orbit Determination Error 

Case IS 

(umber 

Case 1-1 

Case 1-2 

Case 1-3 

GPS positioning 
random error 
[m] 

Along 

0.02 

0.20 

2.00 

Cross 

0.02 

0.20 

2.00 

Radial 

0.04 

0.40 

4.00 

Orbit Determination Error [m] 

8.19 

9.29 

16.32 


The result tells us that the orbit determination error will be more than 8 m even if we use a very accurate 
GPS receiver whose positioning accuracy is only 2.0 cm in the along and cross track, and 4.0 cm in the 
radial. When we use an ordinal GPS receiver such as 2.0 m positioning accuracy the determination error 
increases up to about 16 m. 


5-2-2 Case2: Influence of Initial Coefficients 


To figure out the influence of the initial coefficients to the orbit determination accuracy, three cases with 
different initial coefficient are simulated. Note that the coefficient of the air drag is set to 0.0 because the 
air drag is basically zero except the small part of the orbit near perigee. Table 5-2-2-1 summarises the 
result. 


Table 5-2-2-1: Initial Coefficients and Orbit Determination Error 


Case Number 

Case 2-1 

Case 2-2 

Case 2-3 

Initial 

Coefficients 

Air Drag 

0.0 

0.0 

0.0 

Solar Radiation 
Pressure 

7.06e-10* 

1.77e-9** 

3.53e-9*** 

Orbit Determination Error [m] 

8.19 

31.36 

75.60 


* True value x 2, ** True value x 5, *** True value x 10 


When we see the result we can easily find that the influence of the initial coefficient are very large. If we 
use the value which is ten times larger than the true one the determination error becomes 75.6 m. Perhaps 
the reader of this paper may think that this difference is "too large". However the real coefficients of the air 
drag and the solar radiation pressure are very small so the tens or hundreds times variation is not unusual. 
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5-3 Simulation Result 


5-3-1 Contribution of ACC 

To improve the orbit determination accuracy, an on-board ACC is used. When the ACC is available we can 
get the information of the non-conservative force vectors i.e. the air drag and the solar radiation pressure. 
The information from the ACC allow us to estimate more suitable coefficients, which update the force 
model used in the orbit estimator and improve the determination accuracy. Figure 5-3- 1-1 shows the orbit 
determination error at apogee by the estimation with using the on-board ACC. 


a 

5 1 



Case 1-3 


Case 1-2 
Case 1-1 


3.0e-17 3.0e-14 3.0e-ll 3.0e-8 3.0e-6 3.0e-4 

Accuracy of ACC [m/s2] 

(a) Case 1 


m 

3 



Case 2-3 


3.0e-17 


3.0e-14 


3.0e-6 


3.0e-l 1 3.0e-8 

Accuracy of ACC [m/s2] 

(b) Case 2 

Figure 5-3-1-1: Result of Orbit Determination with ACC 


Case 2-2 


Case 2-1 


3.0e-4 


The horizontal axis shows the accuracy of the on-board ACC and the vertical axis shows the orbit 
determination error at apogee. Figure (a) is the result of Case 1, which shows the errors with three different 
GPS positioning accuracies and the contribution of the on-board ACC to the orbit determination precision. 
Figure (b) is the result of Case 2 showing the errors with three different initial coefficient set. 

At first let's see the result of Figure (a). As you can see the determination accuracy significantly improves 
by the use of the ACC. When we use the hybrid method with an super accurate GPS receiver with 2.0 cm 
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positioning accuracy and a ultra accurate ACC with 3.0e-17 m/s2, then we can determine the orbit in an 
accuracy of 53.0 cm at apogee. Even in the worst case whose determination error is about 16 m when it is 
determined only by GPS, the error decreases down to 5.4 m. 

On the other hand, Figure (b) shows the necessity of the accuracy of the initial coefficient. As it was said in 
the previous section the influence of the initial coefficient to the orbit determination accuracy is very large. 
Even if a very accurate ACC is used for the estimation the determination error still remains more than 30 m 
when the initial value which is ten times larger than the true value is used. In other words, this result tells us 
that we have to determine the initial coefficient as precisely as possible if we need to fit the orbit 
accurately. 

5-3-2 Estimation of the Coefficient of Air Drag 

As you already know, the orbit estimator of this study estimates the two coefficients of the air drag and the 
solar radiation pressure. 



(a) Case 1 


HEO Orbiter: Visibility Analysis 



(b) Case 2 


Figure 5-3-2-1: Result of Estimation of Air Drag Coefficient 

Here the estimation result of the air drag's coefficient is shown. Figures 5-3-2- 1 show the results of the 
estimation with four different ACCs. Similarly to the previous section Figure (a) corresponds to Case 1 and 
Figure (b) is the result of Case 2. The horizontal axis shows the time period starting on 00:00:00.0, 16 th , 
May, 2006. The estimation time period is roughly 4 hours which corresponds to the half of the period of 
ASTRO-G's orbit. The vertical axis shows the value of the estimated coefficient. 

When we scrutinise the results shown in Figure (a) we find that the estimated value is negative number i.e. 
”-9.72e-12" at the end of the estimation time period. Frankly speaking this result is odd because the 
coefficient should always be positive value. 
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However, in this study a very simple atmospheric density model is applied. According to the density model 
the air density has some value when the user satellite's altitude is below 1,750 km, but in the most of the 
part of the orbit this value is zero. In contrast, the orbit estimator simply assumes the coefficient is a 
constant value. Hence the estimator has to make the value negative to cancel the effect of the air drag 
applied to the user satellite near perigee. This is the why. 

Figure (b) also shows a similar result. The three lines appearing in the graph correspond to Case 2-1, Case 
2-2 and Case 2-3 respectively. According to the result all of the three values converge to the same value as 
Figure (a) i.e. "-9.72e-12". On the other hand, the shape of the three lines are different each other while 
those in Figure (a) are perfectly same. This is because of the initial coefficient of the solar radiation 
pressure. The force sensed by the on-board ACC consists of the air drag and the solar radiation pressure, so 
the different coefficient of the solar radiation pressure affects to that of the air drag. 

5-3-3 Estimation of the Coefficient of Solar Radiation Pressure 

Here we see the estimation result of the solar radiation pressure's coefficient. Figures 5-3-3- 1 show the 
result. 

Figure (a) shows us the good estimation result of the coefficient. The value converges to "3.32e-10" while 
the true one is "3.53e-10". The shapes of the lines of those three different cases are completely same, which 
means the influence of the GPS positioning accuracy to the estimation of the coefficients are very small. 

Conversely, Figure (b) shows different characteristic. In Figure (b) the estimated values of the coefficient in 
those three cases are different each other. Obviously this is because of the initial coefficient set. The solar 
radiation pressure's coefficients of Case 2-1, Case 2-2 and Case 2-3 converge to "3.3e-10", "4.5e-10" and 
"6.6e-10" respectively, while the true one is "3.53e-10". The difference of those values are very small, but 
those small variations of the estimated coefficients cause the large orbit determination error as you saw in 
the previous section. 


HEO Orbiter: Visibility Analysis 





(a) Case 1 


Case 1-1 
Case 1-2 
Case 1-3 


HEO Orbiter: Visibility Analysis 



Tie 16 100 200 300 «00 

Hay 2006 Tine 


(b) Case 2 

Figure 5-3-3-1: Result of Estimation of Solar Radiation Pressure Coefficient 
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5-4 Influence by ACC's bias 


In this study, the following ACC bias model is applied as a test. 

Bias = b + b(t ref - 1)+ b(t ref - 1) 1 

At this present time the research for the ACC's bias is not sufficient. However, by using the brief model 
described in the equation (5-4-1) the influence by the bias was tried to figure out. 


Table 5-4-1: Val 

ues of ACC bias 

Bias 

Values 

b 

2.0e-l 1 

b 

2.0e-14 

b 

2.0e-17 


In this study the accuracy of GPS is same as casel-1 in Table 5-2- 1 - 1 . The random noise of ACC is 3.0e-l 1. 
The initial coefficients of the air drag and the solar radiation pressure of the estimator are twice larger than 
the true value. In this case the peak range error became about 3.63m while that of the same case without 
ACC bias is about 0.53m . To reduce this error the ACC bias should also be estimated by the estimator. 

6. Conclusion 

This study revealed the feasibility of the hybrid orbit determination method by GPS and ACC. If only GPS 
measurement data is used for the determination at least 8 m error will emerge at apogee even if a super 
accurate GPS receiver is used. However the ACC allows us to estimate the coefficients of the 
non-conservative forces i.e. the air drag and the solar radiation pressure applied to the user satellite, which 
reduces the determination error down to only 53 cm. This study also tells us that the more accurate ACC 
makes the determination accuracy better, but the contribution of the ACCs which are more accurate than 
"3.0e-8 [m/s2]" saturates. 

The influence of GPS positioning accuracy to the orbit determination error is not small. A few meters 
positioning error at perigee will produce about 16 m error at apogee. This accuracy should be taken into 
account through the use of dual frequencies, a more detailed Kalman filter capable of modelling the clock 
bias model of the GPS receiver when we need to determine the orbit precisely. The use of GPS weak 
signals from the side lobe of the transmission antenna may also be useful to achieve better result. 

The biggest factor which affects to the determination accuracy is the initial coefficients set of the air drag 
and the solar radiation pressure. Generally speaking, these values are so small that the ten or hundred times 
variations of those values are common. However this study revealed the ten times difference will cause 
more than 75 m error at apogee. Hence it is vitally important to set those initial coefficients as precisely as 
possible if a very accurate orbit determination accuracy is required. 

As future works, the problem of an ACC's bias still remains. In real situation the output data from an 
on-board ACC contains some bias which changes over time. This problem should be discussed when this 
method is applied to the application of real spacecrafts. To solve this problem it is intended to estimate the 
ACC's bias as well as other state parameters. The bias model should be considered for the modification. 

The calibration of the on-board ACC is also important problem. This also should be taken into account to 
make this algorithm more practical. 
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